Gray wolf optimizer with bubble-net predation for modeling fluidized catalytic cracking unit main fractionator

Fluidized catalytic cracking unit (FCCU) main fractionator is a complex system with multivariable, nonlinear and uncertainty. Its modeling is a hard nut to crack. Ordinary modeling methods are difficult to estimate its dynamic characteristics accurately. In this work, the gray wolf optimizer with bubble-net predation (GWO_BP) is proposed for solving this complex optimization problem. GWO_BP can effectively balance the detectability and exploitability to find the optimal value faster, and improve the accuracy. The head wolf has the best fitness value in GWO. GWO_BP uses the spiral bubble predation method of whale to replace the surrounding hunting scheme of the head wolf, which enhances the global search ability and speeds up the convergence speed. And Lévy flight is applied to improve the wolf search strategy to update the positions of wolfpack for overcoming the disadvantage of easily falling into local optimum. The experiments of the basic GWO, the particle swarm optimization (PSO) and the GWO_BP are carried out with 12 typical test functions. The experimental results show that GWO_BP has the best optimization accuracy. Then, the GWO_BP is used to solve the parameter estimation problem of FCCU main fractionator model. The simulation results show that the FCCU main fractionator model established by the proposed modeling method can accurately reflect the dynamic characteristics of the real world.

Fluid catalytic cracking is one important part of the petroleum refining processes. The heavy oil is reacted in a lifting tube having a temperature of about 500 °C in the presence of a molecular sieve catalyst. The hydrocarbons of the born can be separated into cracking gas, high octane gasoline and diesel. It is one of the important measures to convert heavy oil into light oil. The catalytic cracking main fractionation column is a master device that achieves product separation during catalytic cracking, which is an important part of catalytic cracking. Therefore, establishing a precise main fractionator model has a great significance to improve oil product quality 1 . The catalytic crack process is full of uncertainty, complexity, nonlinearity and coupling among variables etc. 2,3 . It is very difficult to obtain an accurate mathematical model by traditional modeling methods, they have been unable to meet its requirements of building a high-precision mathematical model 4,5 . The modeling problem of complex system can be transformed into a complex optimization problem. Researchers found that intelligent optimization algorithms can solve the complex optimization problems [6][7][8] . Intelligent optimization algorithms have no strict requirements for optimization problems, and easy to operate and implement 9 .
With the development of artificial intelligence technique, intelligent optimization algorithms are also used in petroleum industry modeling 10 . Li et al. applied Tabu Search to solve chemical industry optimization problems, and obtains the global optimal solution in specific constrained optimization problems by setting different parameters 11 . Tao and Wang proposed an RNA genetic algorithm to estimate the parameters of a heavy oil thermal cracking model and a fluid catalytic cracking unit (FCCU) main fractionator 12 . Li et al. used grey wolf optimizer (GWO) to optimize multi-classification Twin Support Vector Machine to explore concealed reservoirs 13 . Hu et al. used a normal distribution function to improve the search mode of fruit fly optimization algorithm to predict the oil pipeline energy consumption 14 .
Inspired by the preying activity of gray wolf, grey wolf optimizer (GWO) is proposed in 2014 15 . GWO has strong search performance, few parameters and easy implementation. Some researchers used GWO  Nadimi-Shahraki et al. used the dimension learning hunting search strategy to improve GWO, so that each wolf constructs a neighborhood, and the adjacent information is shared among wolves, which alleviated the lack of population diversity, the imbalance between development and exploration and the premature convergence 24 . Liu and Wang proposed a GWO with RNA cross operation to improve the global optimization ability and changed the adaptive parameter to balance the exploration and development ability 25 . The improved methods as above provide the reference values. The global search ability (detectability) and local search ability (exploitability) of GWO needs to study for reaching an ideal balance. The method proposed in this work can efficiently balance the detectability and exploitability. In this work, in order to overcome the shortcomings of basic GWO, we combine the whale search scheme and the Lévy flight for the head wolf α search strategy to propose the novel gray wolf optimizer (GWO_BP), and adopt GWO_BP to estimate the parameters of the FCCU main fractionator model. The highlights of this paper are as follows: • The gray wolf optimization with bubble-net predation (GWO_BP) is proposed.
• In order to enhance the global search ability and accelerate the convergence speed, the bubble-net predation of whale search scheme is applied to update the head wolf α position for GWO_BP. • And the Lévy flight is used to the head wolf α to update the positions of wolfpack for overcoming the disadvantage of easily falling into local optimum. • GWO combined with whale bubble predation and Lévy flight for head wolf α search strategy is applied to the parameter estimation of FCCU main fractionator model. The experimental results show that the proposed modeling method can track the dynamic characteristics well.
The chapter structure is as follows. The second section introduces the basic GWO. The third section proposes a novel GWO (GWO_BP), which merge the whale bubble-net predation search scheme and Lévy flight. The fourth section uses the typical optimization test functions to test the GWO_BP and the basic GWO and PSO. The GWO_BP is applied to estimate the model parameter of FCCU main fractionator to verify the effectiveness and feasibility.

Basic gray wolf optimizer
GWO is a newly intelligence optimization algorithm proposed by Mirjalili et al. 15 . GWO is inspired by the prey hunting activities of gray wolves. Gray wolves are social canine animals that live at the top of the food chain, and they adhere to a strict hierarchy of social dominance. There are four social hierarchies in the gray wolf population α , β , δ , ω . Where α is the leading wolf in the wolf pack, and it makes major decisions about activities in the wolves, In the GWO, α is the best fitness value in the wolves. β wolf obeys the leading wolf and helps him to make decisions. β wolf can dominate the other grades of wolves, and β is the second best solution of fitness value in GWO. δ obeys α and β in the wolf pack, and at the same time dominates the wolves of the remaining level. δ is the third best solution of fitness value. The rest of the candidate solutions are ω , ω wolves usually need to obey wolves at higher social levels. The optimization process of GWO is guided by α , β and δ , after judging the position of the prey which is the optimal solution, α , β and δ lead ω to surround the prey. Finally, the optimal value is found by iterating continuously.
GWO can divide the whole process of hunting prey into three stages, which are encircling, hunting, attacking, and finally capturing the prey, The detailed algorithm is described as follows 15 : (1) Encircling prey: after the wolves lock the location of the prey, they will slowly move to the prey for encirclement. In the encirclement process, the distance between the gray wolf and the prey can be calculated by Eqs. (1) and (2) 15 : where − → D indicates the distance between the wolf and prey, t is the current iteration, − → X P (t) is the location of prey after the iteration t , − → A and − → C are coefficients, and can be calculated as follows 15 : www.nature.com/scientificreports/ where − → r 1 , − → r 2 are random value in [0,1]. − → a is the convergence factor and decreases linearly from 2 to 0 with the iteration proceeds, the definition is as follows 17 : (2) Hunting: the hunting process of gray wolf population is usually guided by the head wolf α , β and δ assisted the wolf α . Therefore, GWO assumes that α , β and δ are related to the possible positions of prey, and updates according to the positions of these three optimal solutions, the expression are as follows 15 : (3) Attack prey: This step is the last stage of the hunting process, the wolves siege and capture the prey (obtain the optimal solution). This stage can be realized by decreasing the value of − → a in Eq. (3). When the value of − → a decreases linearly from 2 to 0, the corresponding value of − → a also changes in the interval [− a, a]. When the random value of A is above [− 1,1], the next position of the wolf may be anywhere between its current position and its prey position, when |A|> 1, wolves are currently moving away from position of prey to find new potential prey. Random parameter C in Eq. (4) ranges from [0,2], parameter C randomly enhances (C ≥ 1) or weakens (C < 1) the influence of the target wolf on the computational distance, this helps to enhance the detectability of the algorithm and avoid local optimum values.
Gray wolf optimizer with bubble-net predation. The bubble-net predation of whales. When the basic GWO encircles the prey, it will gradually make the wolves approach the α , β , δ wolves through continuous iteration, which may reduce the population diversity and lead the algorithm falling into the local optimal values. It is unfavorable for solving the problems with multiple local optimal values. We turn our attention to the whale foraging behavior, and make use of the bubble-net predation of whales to enhance the leading wolf α.
Whales are the largest mammals in the world, and they have a unique method of foraging, called the bubble feeding method. Whales dive into the ocean about 12 m, then create a spiral bubble around their prey, and then swim to the surface. Inspired by this foraging behavior, whale optimization algorithm (WOA) is proposed by Mirjalili et al. 26 . The advantages of WOA are simple operation, fewer parameters to adjust good global optimization ability and the convergence speed is fast. The bubble net predation is expressed by the following formula 26 : represents the length of the ith whale from its prey, that is, the distance between the ith solution and the current optimal solution; − → X * (t) represents the best whale position up to now, − → X (t) represents the current position of the whale; b is a constant, usually b takes 1, which determines the shape of the diameter spiral; l is the random value on [− 1, 1], and '•' represents the dot product.
In addition, the algorithm believes that whales also surround their prey, which is similar to the encirclement mechanism of GWO, the mathematical description is as follows 26 : Suppose we update the whale's position with probabilistic selective contraction encirclement mechanism of P i and probability selection spiral model of Eq. (9). WOA sets when A < 1 , the whale attacked its prey, the mathematical formula is as follows 26 : is the radius of the contraction circle, and the definition of − → A , − → C are as follows 17,26 : When WOA set to A ≥ 1 , the whale is forced to deviate from the prey, so as to find a more suitable prey, which can strengthen the exploration ability of the algorithm and make WOA to conduct global search. The mathematical model is as follows 26 : where − → X rand is a random selection of whale locations (not a pre-optimal solution). In GWO_BP, we will use Eqs. (11)(12)(13) to replace the position update formula of α wolf.
Lévy flight for the head wolf α. In order to improve the population diversity and avoid falling into the local optimal solution, we use Lévy flight to improve the alpha wolf search strategy for global detection.
Lévy flight 27 is a random search method that obeys Lévy distribution and Lévy flight is named after Paul Pierre Lévy, the French mathematician. The search step of Lévy flight is a short-distance and long-distance search alternately. Such a search method has good global search capabilities 28 .
The position update equation of Lévy flight is as follows 27 : where, ⊕ represents point-to-point multiplication, and l > 0 is the step size parameter related to the scope of the optimization problem.
As the equation of Lévy flight is very complex, we use Mantegna algorithm to simulate 29,30 . The calculation equation of step size is as follows: In GWO, the value of α wolf is the closest to the optimal solution, so we add α to s p , a new step size formula is obtained, it is a complex process. The step length calculation formula is as follows: where, µ and ν Obey normal distribution, µ ∼ N(0, σ 2 µ ) , ν ∼ N(0, σ 2 ν ) , σ ν = 1 , The equation of σ µ is as follows 31 : where the value of χ is usually 1.5.
The GWO_BP. From the discussed as above, the flowchart of GWO_BP is shown in Fig. 1, and the procedure of GWO_BP are as follows: Step 1 Set the control parameters of the algorithm. They are the population size S , the maximum number of iterations T max , the dimension of variables to be optimized dim , the sum of upper ub and lower bounds lb in space, and initial population randomly.
Step 2 Calculate the fitness of all individuals in the population, sort by fitness, that is, α = the individual with the best fitness, β = the individual with the second fitness, δ = the individual with the third fitness ranking.
Step 4 Update the position of the wolf pack according to the Lévy flight equation of the head wolf α position in Eqs. (14)(15)(16)(17). Step Step 6 Determine whether the number of iterations reaches the maximum T max or satisfy other algorithm termination conditions. If satisfied, output the current optimal solution, otherwise, return to step 2 and continue.

Experimental results
Test function optimization. In order to verify the effectiveness of GWO_BP, twelve test functions in Tables 1, 2 Table 1 are unimodal functions, they only have one extreme point, which are mainly used to test the accuracy of the local search and the development ability of the algorithm; f 5 − f 8 are multimodal functions, the feature of these functions are that they have multiple local minima, so they also suitable for detecting the global search ability of the algorithm, that is, the detectability of the GWO_BP: whether it has the ability to jump out of the local optimum and search for the  Table 4. Each algorithm of GWO_BP, GWO and PSO is run 30 times for the 12 test functions to calculate the optimal value, mean value and variance. The results of the three algorithms are shown in Table 5, and the optimal values of the three algorithms are marked in bold.    www.nature.com/scientificreports/ Results and discussion. From Table 5, for the comparison of the optimal values of unimodal function f 1 − f 4 , it can be seen that the optimization accuracy of GWO_BP is far beyond several orders of magnitude of GWO and PSO. Therefore, it can be concluded that the local exploration ability of GWO_BP proposed in this work is much higher performance than that of GWO and PSO. At the same time, comparing the mean value and the variance of the three algorithms, the result of GWO_BP is still far better than GWO and PSO. It also can be concluded that GWO_BP is better than GWO and PSO, which has very strong local search ability and good stability; According to the comparison of the optimal values of multimodal function f 5 − f 8 in Table 5, the exploration performance of GWO_BP is better than GWO and PSO. Optimization of GWO_BP for f 6 and f 8 has reached the optimal value of the function, and the average results have reached the optimal value of 0, which proves that the difference between most optimization results and their mean results is small, and the optimization results are relatively stable each time. To sum up, it can be concluded that GWO_BP has strong global search ability and stability; Through the optimization of the function with multi peak fixed dimension f 9 − f 12 , it can be seen from Table 5 that GWO_BP is slightly better than GWO and PSO in the comparison of the optimal value results of the functions.

Modeling of FCCU main fractionator. Description of FCC unit main fractionator.
With the gradual development of national economy, heavy oil catalytic cracking has become a very important problem in today's industrial production 32,33 . The FCC fractionation system is mainly composed of the main fractionator, the overhead oil and gas condensation cooling system, the diesel stripper, the recycle tank and its interrupted reflux. In the typical split flow system process, the high-temperature reaction oil vapor mixture at 450-510 °C comes out from the top of the reactor, entrains a small amount of catalyst powder, enters the desuperheating section of the Table 3. Fixed-dimension multimodal benchmark functions.

Function name Test functions Bounds Dim F min
00030 4028 Table 4. Parameter setting.

GWO_PB GWO PSO
LGWO T max = 500 T max = 500 T max = 500 T max = 500  www.nature.com/scientificreports/ lower section of the main fractionator, contacts with the 250 °C oil slurry countercurrent on the baffle for heat exchange, desuperheates and washes the entrained catalyst powder, and then enters the main body of the fractionator. In the main fractionator, the oil vapor mixture condensed to the saturated state is separated.
Process modeling. We apply GWO_BP to the parameter estimation of the FCCU main fractionator model 34 .
In the main fractionator of a 1.4 million tonnes heavy oil catalytic cracking unit in a refinery, the reaction oil vapor enters the fractionator from the bottom and is cooled and washed from bottom to top. In order to provide sufficient internal reflux, remove the heat in the tower and make the load distribution of the tower uniform, the fractionator is equipped with four thermal cycle systems, namely, top exhaust heat cycle, first middle exhaust heat cycle, second middle exhaust heat cycle and slurry exhaust heat cycle. The factors that affect the dry point are the top temperature, top pressure, top heat discharge and other related parameters, while the factors that affect the pour point are the top load change, the top pressure change, the first and second heat discharge. The change of heat discharge is realized by the change of flow or the change of temperature, so we set the controlled variables are the top temperature y 1 , the crude gasoline dry point y 2 , the pour point of light diesel oil y 3 ; the control variables are the circulating flow u 1 , the first medium flow u 2 , the second middle flow u 3 . 100 input and output data are used for the parameter modeling; The dynamic model of FCCU main fractionation column process is as follows 34 : Due to the coupling between y 1 and y 2 , Zhong et al. established the model of the main fractionation column of the FCCU based on the field data of typical FCC system. Therefore, parameters of FCCU dynamic mathematical model provided by Zhong et al. are used 34 .
GWO_BP, LGWO and GWO are used to estimate the model parameters in the main fractionator and the algorithms parameter setting is shown in Table 4. The modeling results are shown in Figs. 2, 3 and 4.
The model error formula is as follows 35 : where N is the number of samples, ŷ i is the model output, y i is the model actual output. The adjustable range of y 1 ~ y 3 is 0.0-0.8, and the allowable error range is within 0.1. The output error of GWO_BP, LGWO and GWO are shown in the Table 6. From the modeling results in Figs. 2, 3 and 4, it can be seen that the modeling results of this method reflect the actual characteristics of the actual system well, while the modeling results of GWO are far less ideal than that of GWO_BP. From the error analyzed in Table 6, it can be seen that the error of the improved algorithm is closer to 0. In addition, in order to further verify the effectiveness of the new algorithm in model building, the error is compared with LGWO. The experimental results show that the output error of GWO_BP is less than LGWO.

Conclusions
To realize the advanced control of product quality in a refinery is to establish a high precision dynamic model of complex industrial process. In this work, by combining bubble-net predation of whales and Lévy flight, a novel GWO (GWO_BP) is proposed for parameter estimation of FCCU main fractionator model. The new algorithm can make up for the imbalance between exploration and development of the original GWO. On the one hand, the whale bubble predation method is replaced by the surrounding predation method of the head wolf α in GWO to enhance the global search ability; On the other hand, the Lévy flight improved the head wolf α search strategy is used to iteratively update the wolf swarm to overcome the disadvantage of the algorithm falling into local optimization, so as to speed up the convergence speed and improve the convergence accuracy of the algorithm as a whole. GWO_BP, GWO and PSO are applied to 12 typical test functions. The experimental results show that (18)  www.nature.com/scientificreports/ the performance of GWO_BP is much better than the others. And compared with basic GWO and LGWO, the results reveal that the FCCU main fractionator model predictive outputs of GWO_BP are in better agreement with the actual experimental data.